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ABSTRACT 


A 

A  least  squares  estimator  4  is  derived  for  the  state  transition 
matrix  4  of  a  linear,  stationary  sampled  data  system  operating  in  a 

A 

stochastic  environment.  The  estimator  4  is  shown  to  be  unbiased  and 
minimum  variance  under  the  condition  of  full  observability  of  the  state 
vector  of  the  system.  The  estimator  is  also  shown  to  be  the  Maximum 
Likelihood  Estimator  for  the  case  of  the  stochastic  environment  having 
Gaussian  statistics .  The  estimation  scheme  is  compared  with  two 
other  recently  published  estimation  schemes ,  both  of  which  are  shown 
to  be  special  cases  of  the  scheme  herein  presented. 
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1 .  INTRODUCTION 


The  problem  to  be  considered  is  the  identification ,  in  a  stochastic 
environment,  of  linear,  stationary  (constant  coefficient),  sampled  data 
process  dynamics.  The  problem  is  more  concisely  defined  as  follows: 

Given: 

1)  A  system  whose  behavior  may  be  described  by  a  set  of  linear, 
constant  coefficient ,  differential  equations . 

2)  A  statistical  description  of  the  excitation  of  the  system  (the 
vector  w  in  the  figure  below) . 

3)  A  sequence  of  observations  on  the  state  vector  (x^  in  the 
figure) . 

Problem: 

Determine  a  "best"  estimate  of  the  system  dynamics. 


SYSTEM  WITH 

A 

x(t)  j 

w 

UNKNOWN 

DYNAMICS 

A. 


Figure  1 .  Schematic  representation  of  the  problem 

Identification,  defined  in  this  manner,  amounts  to  estimating  the 
solution  of  the  differential  equation  postulated  in  1)  above. 

This  problem  may  appear  at  the  outset  to  be  a  rather  restricted 
one.  It  becomes  apparent,  however,  that  a  study  of  the  problem  takes 
on  added  significance  when  viewed  in  the  context  of  a  more  general 
question.  Given  a  vector  of  random  time  functions  produced  by  some 
unknown  process,  can  the  process  be  adequately  represented  by  a 
linear,  stationary  model,  and  if  so,  which  linear  model  is  the  best  rep¬ 
resentation  of  the  system?  This  work  is  motivated  primarily  with  the 
hope  of  finding  application  in  the  investigation  of  such  questions . 

Significant  contributions  in  this  area  have  been  made  by  Ho  and 
"Whalen  [23  and  Lee  [33.  The  purpose  of  this  work  will  be  to  extend 
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the  methods  suggested  by  the  above  authors .  The  approach  here  wilt  be 
to  derive  a  scheme  which  minimizes  a  quadratic  cost  criteria,  subject 
to  some  rather  limiting  assumptions.  An  attempt  will  be  made  to  analyze 
the  consequence  of  relaxing  these  restrictive  assumptions.  Finally, 
the  methods  of  Ho  and  Whalen  [2  ]  and  Lee  [  3  ]  will  be  generated 
from  the  method  herein  derived.  Acknowledgement  Is  given  to  LT  Ralph 
E.  Hudson,  USN,  who  first  suggested  the  algorithm  to  be  presented. 

2 .  BASIC  MATHEMATICAL  MODEL 

rA 7,  . ,  ,  .  s  .  , 

Consider  a  dynamic  system  whose  behavior  may  be  defined  by  a 
set  of  n  linear,  constant  coefficient,  differential  equations. 

x  =  F  x  +  Dw* 

The  solution  is  given  by 

t 

x(t)  =  4  (t-t  )  x  +  f  4(t-r)  Dw*(t)  d  T 

O  O  “ 

o 

^  if  i  va  ./  ;■■■■ 

where  $  (t-t  )  satisfies 
o 

£  *  “-‘o’  “  F  4  “-‘o’ 

,  .  ■  O  ,  '  ■  !>'  H.  < 

Introducing  a  sampling  device  of  period  T  and  a  zero  order  hold 
on  the  excitation  signal  w*(t)  makes  possible  the  representation  of  the 
system  at  multiples  of  the  sampling  instant  T  by  ' 

Vi  ■  4xk  +  rV 

where 

x^  ■  x(kT) 

4  ■  4(T) 

f3  i  ni  -  •/;  -  *  *’■:  T  »\ 

r  TW-J  4(T-T)DdT, 

O 

utIJS  Oi'i  1<'  ,'u  b ''T.r-  di.Sj  C.  .t'l.Hi'w  ‘J  O' 

i.Mte.tx'a  o3  *  -l  ii  w  jhow  uriJ  i<'  ,  L  j  :>J  > 
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n  ;j 


■■id , 


.-m. 


,  M 


Figure  2.  Schematic  representation  of  the  difference  equation 


3 .  DERIVATION  OF  THE  IDENTIFICATION  ALGORITHM 
Consider  a  process  of  the  form 


Vi 


4  x,  +  w, 


(3.1) 


wk  'rwk 

Here  is  defined  to  be  a  gaussian  sequence  of  zero  mean,  with  co- 
variance  matrix  Q.  Suppose  that  a  sequence  of  m+1  observations  are 
made  on  the  state  vector  x;  (xq,  ,  ,  ,  .  x^) ,  from  which  it  is  desired 
to  form  an  estimate  of  the  state  transition  matrix  4 ,  assumed  to  be  un- 

A 

known.  The  so-called  "  lea  sb  squares"  estimator  is  that  estimator  4 
which  minimizes  the  quadratic  cost  function  L  (a  scalar),  where 


L  ■  E .  <Vi  -  *  xk)T(Vi  -  *  xk) 

k*“0 


(3.2) 


Equivalently,  4  minimizes  the  trace  of  J,  where 
n-lho  *  at 

Jm? \  <xk+i  -  *  V'Vi  - 4  xk> 

k«0 


(3.3) 


From  inspection,  L  *  tr(J). 

>  ?  '  i  Ofi"  St'S’  ..  1  v  r).  ,  v  : 

The  classical  method  of  solving  a  problem  of  this  nature  is  to  set 

A 

the  gradient  of  the  cost  function  with  respect  to  the  estimator  4  equal 

v  n  •  :  -  v)  \  >  ■ 

I  'Li 


to  zero.  The  solution  to  this  gradient  equation  is  sufficient  to  es¬ 
tablish  a  minimum  cost.  The  method  cannot  be  applied  here,  however. 

A 

since  the  gradient  operation  is  not  defined  for  the  matrix  $  . 

To  find  a  solution,  it  is  proposed  to  reformulate  the  minimization 
criteria,  and  then  show  the  subsequent  solution  also  minimizes  the 
original  cost  function,  tr(J). 

Let  the  order  of  the  system  be  n,  and  partition  the  state  transition 

T 

matrix  4  into  n  row  vectors  ,  i*l,  2,  .  .  .  n. 

Then  the  process  model 

Vl  “  4xk  +  wk  (3-4) 

may  be  written  as  m  scalar  equations , 

Vi  “  *1  xk  +  wk  •  (3-5) 

Here  the  superscript  denotes  the  scalar  component  (or  mode) 

of  the  vector.  From  the  m+1  vector  observations,  (x  ,  x, ,  .  .  .  x  ), 

o  l  m 

define 


•  .! 

xT  -  [Vl.  .  .^.j] 

(nxm  matrix) 

(mxl  vector)  (mxl  vector) 


From  (3.5),  these  newly  defined  terms  are  related  by 


x  <t>i  +  ct 


(3.6) 


a  . ;  -  i  .  , 

Now  consider  an  estimator  ,  which  minimizes  the  quadratic  cost 

1  /».  ,  \  be rfjeA:  it'  •  i 

function  g.(a  scalar),  where 

#  X'aM  A  T  -  ’»AU  ?<::j  *}■■'  2*1,  *♦ 

gl«(v1-X0i)A(v1-X^i)  (3.7) 
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Comparing  equations  3.3  and  3.7 ,  it  is  seen  that 


9i  "  jii  '  where  1  “  <jik* 


Thus,  the  sum  of  the  new  cost  function  g^,  i*l ,  2 ,  .  .  .  n,  is  pre¬ 
cisely  the  same  as  that  previously  defined,  in  that 
n 

I  g.  -  tr<J)  (3.8) 

i-1 


A  sufficient  condition  for  the  minimisation  of  tr(J),  then,  is  that  g  be 

A  * 

minimised  for  i*l,  2,  .  .  .  n.  Further,  g  is  not  a  function  of  0  , 

.  1  A  * 

if) .  This  implies  that  the  optimum  estimator  0  does  not  depend  upon 

- :  i  '  f  "  ■  ‘  4 

the  choice  of  the  estimator  selected  for  the  jth  row  of  the  matrix  es- 

A  A 

timator,  4  .  Thus,  each  ^  may  be  found  independently. 

To  accomplish  the  minimisation,  expand  equation  3.7  and  perform 
the  gradient  operation. 

7  .  »!  “  -<XTVj)T  -  vj  X  +  2  i*  XTX  (3.8) 

*1 

A 

Setting  VA  g1  *  0  and  solving  for  ^ 


itT  -  v*  X(XTxf 1 


(3.9) 


T 

Equation  (3.9)  assumes  X  X  to  be  nonsingular. 


The  conditions  for 


nonsingularity  are  discussed  in  Section  5. 

A 

The  optimum  estimator  4  can  now  be  formed  by  placing  the  row 

A  f 

vectors  0  ^  in  matrix  form . 


•A 

4 


X(XXX)_1 


(3.10) 


Finally,  writing  4  as  a  function  of  the  original  sequence  of  vector 

observations  x  ,  x, ,  .  .  .  x  ,  gives 
o  l  m 


»  n>-l  m-1 

*  ‘  E  Vi  \  E  Vk 

k*0  *  1  K  k*0  *  * 


-1  r:  ,, 


(3.11) 


Here  £  x,  x,  is  recognized  as  a  finite  approximation  to  the 
k-0 

autocorrelation  function,  R(kT),  evaluated  at  k«=0.  Similarly, 
m-1  T 

L  xjc+i\  is  a  approximation  to  R(T).  Making  this  identification 
k“0 


Ri  ®0> 


(3.12) 


where 


Rj  *=  R(T) 


Ro  “  R(o)  • 


4.  DEVELOPMENT  OF  A  RECURSIVE  FORM 

A 

Note  that  (3.12)  gives  the  optimum  estimator  4  ,  based  on  a 

sample  of  size  m+1.  The  estimator  is,  then,  a  function  of  m,  allowing 

'  ^ 

the  functional  notation 


!0*  ctr  >>;  fDAc. 


I  ■  R  (m)  [R  (m)T1 
ml  o 


To  keep  the  index  consistent  with  the  sample  size,  define 


i^O'l  (jiij 


R  (m)  ■  2  x  .  .x 

X  k,0  Kfi  x 


(4.1) 


R  (m)  ■  £  xx 
°  k-0  *  K 

Incrementing  the  index  in  (4.1)  gives 


^m+1  ■  CRjfrn+1)]  L  R0(m+1)] 
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From  the  definition  of  Rj  (m) 


Cl  “  CRlfm,  +  xm+lxmKRo(m)+xmxmrl  (4'2> 


To  ease  the  reader  through  the  remainder  of  the  recursive  development. 


x  Px  x  »x  , .  R  *  R  (m) 

m  1  m+1  o  o 


Then 


Ci  ■  <vx1xT)(r0+xxT)'1 


Appealing  to  the  matrix  Inversion  lemma,  [  5  ] 


(R  +xxTf 1  -  R- 1  (I  -  xCxT  R_1  x  +  l]”1  xT  R"1) 
o  o  o  o 


(4.3) 


Recognizing  R,R_1  as  ^ 

1  o  m 


4m+i“ci+xix  roki-x(x  r;  *+i>'  *  r;i 


1  -l  -l 

Define  the  scalar  a  =  (x  R^  x_  +  1) 

m  mom 


Collecting  terms  and  simplifying 


i  ,,  -  +  o  (x  -  *  x  )xj  r"1^) 

m+l  m  m  m+1  m  m  m  o 


(4.4) 


Define  P  *=  R  (m) . 
m  o 

From  equations  (4.3)  and  (4.4),  the  final  recursive  forms  can  be 
written  down. 


4  *4  +  a  (x  x  )xT  P 

m+1  m  m  m+1  m  m  m  m 


(4.5) 


o  r^o.  t  v. •  y  , 


n  £  m 
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(4.6) 


m+1 


P  (I  -  a  x  x  p  ) 
m  m  m  m  m 


a  ={x  P  xT  +  i)  1 
m  m  m  m 


(4.7) 


Implicit  in  the  use  of  these  recursive  forms  is  the  utilization  of 

(4.1)  to  generate  an  initial  estimator  £  and  an  initial  matrix  P  .  This 

m 

fact  would  seriously  complicate  a  mechanization  of  the  iterative  process, 

However,  this  difficulty  can  be  avoided  by  the  appropriate  assignment 

of  values  to  the  initial  estimator  4  ,  and  the  initial  P  matrix.  Digi- 

m 

tal  simulation  has  produced  no  apparent  degradation  of  the  estimation 
scheme  when  using  this  ruse.  This  subject  is  further  discussed  in 
Section  10. 

5 .  NECESSARY  AND  SUFFICIENT  CONDITIONS  FOR  IDENTIFIABILITY 
Equation  (3.9)  assumed  the  nonsingularity  of 
m-1  „ 


l  =  E  x,  x 
o  ,  n  K  K 

k=0 


(nxn)  matrix 


Consider  the  decomposition  of  Rq  =  X  X 


o 

2 


n 


x 


.  x 


m-1 


n 

cm-l 


(nxm)  matrix 


Here  the  superscript  again  denotes  the  element  or  mode  of  thq  state 
vector. 

From  the  fact  that  X  must  be  of  rank  n  If  R  is  to  be  non- singular, 

o 

it  is  immediately  apparent  that  m  in  is  a  necessary  condition  for 
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identlflabllity .  The  number  of  observations  on  the  state  vector 

must  be  equal  to  or  greater  than  the  system  order,  n.  From  Inspection 

of  X,  It  is  seen  that  if  any  mode  of  the  system,  x^  ,  k  =  0,  1 ,  .  .  ,  m-1, 

remains  at  zero  for  all  k,  or  if  any  two  modes  remain  at  some  constant 

value,  X  will  be  of  rank  r  <n,  and  R  will  be  singular.  These  con- 

0  -1 

ditions  imply  that  all  modes  of  the  system  must  be  excited  if  Rq  Is 
to  exist. 

Finally,  to  insure  sufficiency  of  tie  above  conditions,  recall  that 
the  solution  of  an  n  order  differential  equation  generates  exactly  n 
linearly  independent  solutions .  The  solution  for  the  n  modes  of  the 
system,  then,  form  sets  of  linearly  independent  vectors ,  provided  all 
modes  are  present. 

T 

From  the  fact  that  Rq  has  the  decomposition  X  X,  it  follows  that 

it  is  positive  semidefinite .  If  the  system  is  identifiable,  R  will  be 

o 

nonsingular  and  therefore  positive  definite.  R  is  further  seen  to  be 

T  0 

symmetric,  from  the  decomposition  X  X. 

6.  IDENTIFICATION  OF  A  SYSTEM  WITH  CONSTRAINED  OBSERVABILITY 
The  material  presented  in  this  section  is  an  extension  of  a  de¬ 
velopment  of  Lee  [  3  ]  ,  who  considers  the  case  of  a  scalar  observable. 
Consider  first  a  free  (unforced)  linear  process  of  order  n. 

xk+l-txk"*k+lxo  (6-1) 

Let  the  observability  of  the  state  vector  be  constrained 

*k  -  Hxk  (6.2) 

Here  z^,  the  observation,  is  an  (4x.l)  vector,  l<n,  and  H,  the 

observability  transformation  is  of  dimension  (ixn)  ,  Such  a  system 

is  said  to  be  observable  if  the  initial  state  vector  x  may  be  deter- 

o 

mined  from  the  sequence  of  observations  z  ,  z, ,  .  .  .  assuming 

o  i 

that  H  and  ♦  are  known.  Forming  this  sequence  of  observations  into 

i 

an  (nxl)  vector  gives 
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J  i 


(6.3) 


i 


where  z  denotes  that  observation  required  to  make  y  of  dimension 
v  o 

(nxl) ,  For  example ,  if  n«5  and  l  **2 


Here  z  denotes  the  element  z.  . 

V  z  i 

From  (6.1)  and  (6.2) 

•  l  \  M  ..  1  1  )  '»•"  ’ 


z 

'  ■  d 

H  x 

~H  o’ 

0 

0 

I .  ) 

z. 

H  4  x 

H  4 

1 

o 

J 

• 

• 

• 

y_  * 

*  jr 

m 

m 

1*  * 

o 

• 

• 

• 

•s  j) 

\ 

• 

• 

• 

i  U  .  H. 

z 

H  4  ^  x 

i f  Pt 

H  4^ 

_  V_ 

_  V  o 

J  l>  _ 

x  ■  A  X 
o  o 


V 


tfo  lo  j  POUCIPV.  T6. 


■  I. 


.  * 


;  isb  i'-: 


(6.4) 


H  *  is  defined  in  the  same  context  as  z^  above,  and  may  be  written 
out  for  the  example  of  m*5  >  l  «2 . 

j.  o 


l&ado  Jo 


O-Jiii 


k  osot 

2 


H  4* 


H„« 


!J)d3  aa 


iV.'.otxi  &  bftft  H  nsrr 


•  •  ’>  v  I  •  ioJ  • 


B 


Referring  to  (6.4)/  it  is  seen  that  the  unique  determination  of  the 

vector  x  from  the  sequence  of  observations  z  ,  z, ,  z  ,  is  de- 

o  o  1  v 

pendent  upon  the  nonsingularlty  of  the  matrix  A.  If  A  is  nonsingular, 
H,  4  are  said  to  constitute  an  observable  pair. 

Assume  H,  4  are  an  observable  pair,  and  consider  the  forced 


system 


Vi“*xk  +  rwk 


Define  a  new  vector 

yk‘Axk 

where  the  transformation  A  is  defined  by  (6.4). 
Incrementing  the  index  gives 

yk+i=Axk+i=A(*  xk+rwk> 


(6.5) 


(6.6) 


(6.7) 


Since  A  is  by  assumption  nonsingular 

yk+i“A*  A'lyk+Ar  wk 


(6.8) 


Define  4  *  ®  A  4  A" 


(6.9) 


From  the  definition  of  (6.9),  4  and  4  *  are  similar  matrices, 
and  will  therefore  have  identical  eigenvalues.  Further,  if  the  process 
is  constrained  with  an  observability  matrix  H 


HxJc*HA"1  yk  *H*  yy 


(6.10) 


From  (6.10),  note  the  sequence  of  observations  z  ,  z.  ,  .  .  . 

St  O  1 

will  be  identical  if  generated  by  either  of  the  processes 


yk+i-**yk+Ar  wk 


*k  “H*yk 


(6.11) 


*k+l"*  Xk+I'  Wk 


Zk  “Hxk 


(6.12) 
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4  *  and  H*  can  be  shown  to  be  of  special  form.  To  demonstrate , 
consider  again  a  5tn  order  system  with  two  system  modes  observable. 
That  is,  ^  a  (2x1)  vector,  and  a  (5x1)  vector.  From  (6.9), 


♦  *  *  A  ♦  A 


H  *V  H*V 


where,  as  before, 


4  H  4 


H  *V~ 
V 

_  V  ~ 


(6.13) 


Taking  4  Into  the  left  bracket  of  (6.13)  gives 


4*  -  H  4 


H  4  ■ 


where 


[H-1  C-!3 


(6.14) 


(6.15) 


(  j) 
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From  equation  (6.15) 


C  H_!  H] 


1  0 
0  1 


Lc_i  c] 


1  0  0 
0  1  0 
0  0  1 


•  f :  i 


f . 


[H_1  C3 


0  0 
0  0 
0  0 


Cc_l  H] 


0  0  0 
0  0  0 


performing  the  multiplication  in  (6.14)  gives 

0  0  | 


** 


0  0  j  I 

0  0  I 


E  H-i!  E  C-1 


(6.16) 


In  general,  by  straightforward  extension  of  the  above  procedure 


(•.  ■ 


f  ir  b 


v  'K-Jr./v 


f! .  I 


0  .  .  .1 


T 

*£*+ 1 


T 


(6.17) 


H*  can  be  found  using  (6.10)  and  (6.15),  for  the  example  case 
of  H  a  (2x5)  matrix. 


H*  *  H  A-1  -  H  [  H_x  «  C_1]  - 


5  Olir  9B  r'  '  '  i.t’ISi  it 


oeriUlu 


u  i 


1  0  0  0  0 

0  10  0  0 

V  i  <m 


i i< 
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For  the  general  case 


r  10...” 

1  !  • 

_  l  * 


(6.18) 


where  H*  48  of  dimension  (lx n). 

Reconsider  equation  (6.11)  for  the  special  case  where  the  first 


(n-jt)  elements  of  the  vector  ATw*  are  identically  zero  for  all  w*. 


y.. ,  “4*yu+Arw* 
k+1  k  k 


■k  -H*yk 


A  moment's  reflection  of  the  special  form  of  4*  and  H*  will  show 


(6.19) 


Thus,  for  this  special  case  (the  conditions  being  on  AT),  the  com¬ 


plete  state  vector  of  an  equivalent  linear  system  ♦  *,  H*,  may  be 


constructed  from  the  observations  z^,  provided  only  that  H,  ♦  form  an 


observable  pair.  The  case  of  nonzero  elements  in  the  vector  ATw* 


is  discussed  in  Section  8. 


7 .  SOME  PROPERTIES  OF  THE  ESTIMATOR 


This  section  is  concerned  with  some  properties  of  the  estimator, 
4  .  Recall  (3.6)  and  (3.9), 


vi‘x*i  +  ‘t 


(3.6) 


A  T  T  -1  T  -1 
\  ■  V*  X  (XAX)  -  vA  X  RoA  . 


•v 


(3.9) 


r  (' 

This  formulation  fits  nicely  into  the  framework  of  classical  re- 

j  ■  r  >  < ,  1 1  1 " 

gresslon  analysis  [4  ]  .  This  fact  can  be  utilized  to  determine  several 


■'X 


A 

statistical  properties  of  the  estimator,  -0  .  In  particular,  the  ex¬ 
pected  value  of  the  estimator  may  be  calculated.  Recall  the  vector 
was  formed  from  a  scalar  sequence  of  zero  mean. 


E(wJ)  =  0 


From  (3 . 6)  and  (3 . 9) 


E(vt)  *  E(X  *  X  0i 


(7.1) 


E(  dtT)-  E(v*  X  (X'X)’1)  -  <t>J  (7.2) 

If  has  non  zero  mean,  say  E(  s  6  ^  then 

E(  0{)  -  01  +  X  R"1  (7.3) 

A 

A  requirement,  then,  that  0  be  an  unbiased  estimator  is  that  the 
excitation  be  of  zero  mean. 

To  calculate  the  covariance  of  the  estimator 

.  :.c  - ..  .  J  •.  .  ->  i  .  '  '  ...  . 

cov(  it)  -  E  (  it  -  E  0t)(  -  E  0±)T  (7.4) 

note  that  1 

it  -  E  -  (XTX)“1XT(vi  -  E  Vl)  (7.5) 


v 


i 


cov(vt) 


E(«^T) 


'i*1 


(  O 


2 


a  scalar) 


(7.6) 


airfi  .  (*>' ,  fi  p.-  j<u>  t>  < 

If  is  formed  from  the  sampling  of  a  band  limited  spectrum, 
(7.6)  will  not  hold,  since  it  assumes 


E(wkw|)  ’ 


J-k 


X' 
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2 

a  is  further  seen  to  be  equal  to  q^,  where  Q  *  covfw^). 

Assuming  wj^to  be  from  a  white  spectrum  gives,  using  (7.5)  and  (7.6) 

cov(  0  )  =  (X  R"1)1  E(  C.  C.T)  X  R"1  -  o  2  R"1  (7.7) 

1  o  i  i  o  l  o 

For  the  case  of  a  gaussian  sample  of  zero  mean,  one  can  write 
out  the  joint  distribution  (likelihood)  function  of  the  sample.  Desig¬ 
nate  the  log  of  this  function  by  L(  0^  ,  a2). 

L(  0A  ,  a2)  «=  -§m(log  2  +  log  a*) 

-  i(Vj  -  X01)T(vi  -  X01)/ff12  (7.8) 

The  maximum  likelihood  estimator  is  formed  by  setting 

\L(  V"? ,“° 

Note  however,  from  equation  (3.7),  that  the  gradient  of  the  original 

2 

cost  function  g. ,  is  precisely  the  gradient  of  L(  0.  ,  o  ),  Clearly, 

A  *  1  * 

then,  04  is  the  maximum  likelihood  estimator  in  the  case  of  gaussian, 

it.'.  ■  i 

zero  mean  excitation. 

To  find  the  maximum  likelihood  estimator  for  the  variance  of 
2 

the  excitation,  o,  ,  set 
i 

*  a  L(  V'?>-° 

,  t  «i2-<v1-X*A1-X»1)/»  ,  : 

Replacing  0,  with  0,  in  (7.8),  and  using  (3,9),  this  expression 

1  1  1  ,  n  '• 

reduces  to 

•  i  f  i 

A  2  ,  A  V.T  , 

a  i  *  (Vj  -  0 1  X)  vt/m 
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(7.9) 


This  same  estimator  can  be  derived  from  the  sample  variance  of 
w^,  since 

A  m_1  T  m”l  T 

Q  •  =.0  (wk  wk)/m  “£0  (xk+1  -  txk,(xk+l  -  4xk> /m  (7-10> 

A 

Replacing  4  with  4  ,  and  using  the  fact  that 

4xTx4T  -RliT 

gives 

a  m”1  T 

Q  «  £  ^(w.  wM  «  (R  -  4  R,  )/m  (7.11) 

k=0  °  L 

This  estimator  must  converge  to  Q,  since  4  is  unbiased.  Note 
that  the  same  result  for  (q^)  is  given  by  (7.9).  For  the  case  of  T  an 
n  vector  and  w*  a  scalar,  (7.11)  can  be  used  to  estimate  r  ,  for  in 
this  case, 

q  «E(rw*  w£rT)  «o 2  rrT  (7.12) 

The  properties  of  the  estimator  4  are  summarized  as  follows . 

A 

1)  4  is  the  minimum  variance,  unbiased  estimator  of  4  ,  provided 
the  mean  of  the  excitation  Is  zero. 

2)  For  the  case  of  gaussian,  zero  mean  excitation,  4  is  the 
maximum  likelihood  estimator  of  4  . 

3)  For  the  case  of  gaussian,  zero  mean  excitation, 

‘  >U  ‘Vi;  \  v<:>  »'/t  5  .. 

cov  (  0t)  -  Oj  R”1 

Q  M  (Ro  -  4  Rj  )/m 

8 .  IDENTIFICATION  WITH  NOISY  AND  CONSTRAINED  OBSERVATIONS 
Consider  first  a  system  whose  complete  state  vector  may  be  ob- 

0  /?  *’  A  "**  »  - 

served,  but  the  observations  are  contaminated  by  additive  noise.  Let 
the  observation  noise  be  from  a  gaussian,  zero  mean  distribution. 
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Vi 


=  tx.  +  w. 

k  k 


*k  ■xk+rk 

In  consonance  with  the  definitions  of  suction  3,  define 
i 

ri 
i 


,  i 


n 


(ro  rl  •  •  •  Vl> 

(nxm  matrix) 


(mxl  vector) 

Defining  v^,  X,  and  as  in  section  3  gives 


Vj  -  X  +  ct 


Introducing  the  new  variables 


“i  i 


V.  +p 


W  -X  +  R 


gives,  from  equation  (3.6) 


(8.1) 


(3.6) 


ut  -  o  t  18  (W  -  R)  0 1  +  c  t 

f ,  ■  £  U 


*  .  ,  f 


(8.2) 


Note  here  that  the  arrays  u^  and  W  may  be  formed  from  the  sequence 
of  noisy  observations ,  z^. 

Rearranging  (8.2)  gives 

Vw#l  +  c1  +  p1-R01 
Now  define  a  new  vector 


•do 


il't  10.1 


’9. 


1  1  1  .1 


'  \J  l  t'U 


A  typical  element  of  this  vector  is  given  by 
*k  i  i  T 
{i  -  wk+rk-rk-i«i 


>n  ito 
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From  the  assumed  distribution  of  the  random  vectors  r  and  w 
E<rkr*  )»0  for  j  ^  k 

t 

E(w^  tj)  *  0  forallj,k 

Thus  ^  is  the  sum  of  independent,  zero  mean  gaussian  vectors, 
and  is  therefore  gaussian  with  zero  mean.  The  variance  of  the  vector 
is  the  sum  of  the  variances  of  its  component  parts , 

«*.«?>— i,1-  1 

and  from  (8.2) 

u1-W01T+4i  (8.3) 

From  the  form  of  equation  (8.3),  and  the  conditions  on  £  it 
is  apparent  that  ail  the  proceeding  analysis  leading  to  an  optimum  iden¬ 
tification  scheme  is  also  valid  for  the  case  of  noisy  observations.  The 
deleterious  effect  of  the  measurement  noise  is  readily  apparent  in  the 
increased  covariance  of  the  estimator. 

cov(i)-c*  (V^W)”1 
1  *i 

The  formulation  also  applies  to  estimation  with  constrained  ob- 
servatlons,  as  discussed  in  section  6. 

Vi”*‘yk+Arwk 

■k  -H*yk 


Here  z^  is  an  (  l  xl)  vector,  i< n.  Making  the  assumption  that 
the  first  (n-4)  elements  of  the  vector  A  Tw*  were  identically  zero, 
it  was  shown 

Tk  1 

xk+l 


(  .*) 


(6.19) 
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For  the  more  general  case  where  A  T  w*  has  non  zero  elements , 

th  * 

consider  again  the  example  of  a  5  n  order  system  with  2  observables. 

Let  r  be  a  vector  distribution  matrix,  and  w*  a  scalar.  Let  the 

1  2  5 

elements  of  the  vector  A  T  be  designated  d  ,  d  ,  .  .  .  d  .  From 
equations  (6.11)  and  the  designation  of  A  T  ,  and  utilizing  the  form  of 
H*,  **  , 


zk  + 
2 


d*w* 

k 

d^w* 

k 


o 

d  w.* 
k+1 


“k 

Zk+1  +  d^w* 
k 


k+1 

1 

6k+2 


Define 


- 1 

N 

**  •- 

_ 1 

d*w*  +  d?w*  . 
k  s~  v  .  k+1 

2 

Zk 

At  ‘ 

k 

k+1 

rk“’ 

d^w* 

k 

2 

A 

Zk+1 

u 

1 

Zk+2 

0 

(8.5) 


i .  t  ■ 


Thus 


1 

'k+1 


i  s  >>  > 


*  *  yk  +  a  r  w.  K  , , 


.’.U  I 

(8.6) 


P.  ■  Yu  +  r* 
k  k  k 

From  the  definition  of  r*  it  is  seen  that  the  vector  will  have 

k 


correlation  products  of  the  form 

E(  r*  r*^  )  /  0  for  J  -  1,  2,  ...  n 


(8.7) 
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The  spectral  density  of  r*  is  clearly  a  function  of  A  T  .  Comparing 
(8.6)  and  (8.1) 

(8.6) 

Vi  ■  4xk  +  wk 

(8.1) 

Zk  =Xk+rk 

it  is  noted  that  the  format  of  the  two  sets  of  equations  is  identical. 
Thus,  the  net  effect  of  forming  "pseudo"  state  vectors  from  the  se¬ 
quential  observations  of  dimensionality  l<n  is  seen,  from  (8.6),  to 
be  equivalent  to  adding  correlated  (colored)  measurement  noise  r*  to 
the  state  vector  y^  of  an  equivalent  linear  system,  4*  *  . 

The  correlation  products  of  r*  produce  an  apparent  impasse  in  the 
application  of  the  identification  methods  herein  presented .  The  prop- 

A 

erties  of  0^  as  developed  in  section  7  are  dependent  upon  the  assumption 
that  the  perturbing  sequence  have  the  property  of  statistical  indepen- 

A 

dence.  The  unbiased  and  minimum  variance  properties  of  then,  do 
not  necessarily  hold  in  the  case  of  constrained  observability.  Com¬ 
putational  investigation  has  verified  that,  in  general,  ^  *  will  not  be 
an  unbiased  estimator. 

Arriving  at  this  same  juncture,  but  through  a  different  approach, 

Lee  [  3  3  suggested  using  only  every  n^1  observation  in  forming  the 
estimator  4  *  .  However,  from  (8.5),  it  is  apparent  that  no  amount  of 
time  separation  in  the  observations  utilized  will  render  r*  an  un- 
correlated  sequence.  Computational  experimentation  has  verified  this 
assertion,  and  is  demonstrated  in  section  16. 

Another  approach  considered  was  to  extend  the  order  of  the  es¬ 
timator  to  include  a  coloring  filter.  Since  the  measurement  noise  is 


rk+l 


*  *  y, 


+  a  r  w* 
k 


h  +  rk 
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lumped  with  the  excitation  to  form  4^  in  (8.3),  and  since  there  must 
exist  a  linear  filter  or  order  (n-1)  that  will  produce  r*  from  a  white 

A  * 

spectrum,  it  was  reasoned  that  an  estimator  4**  of  order  2n-l  might 
have  the  necessary  combined  characteristics  of  the  system  being  iden¬ 
tified  and  the  coloring  filter.  This  line  of  reasoning  also  appears  to  be 
invalid,  however,  as  is  demonstrated  in  section  10. 

9 .  A  COMPARISON  WITH  OTHER  ESTIMATION  SCHEMES 

Ho  and  Walen  [  2  ]  have  presented  the  recursion  formula 

fn+1  _  +  (1/X  m,(xm+l ' xm>  xm  (9'1) 


where  X 
vergence 


is  a  scalar  constant.  Using  this  formulation,  the  con- 


lim 

m**-8 


4 


■f ' 


can  be  demonstrated  with  probability  1,  provided  an  appropriate  value 
of  X  is  .selected . 

'  k  ■*  1  - ! '  * 1  5  .  •'  -  '  J  ■  1  '  /  ■  i  ♦  ;  ‘  i 

Compare  this  formulation  with  (4.5) 
a  a  a  t 


4*4  + «  (x  .  -  ^  x  )x*  P 
m+1  m  m  m+1  m  m  m  m 


(4.5) 


on  .  , )  " 


a  ■  (x  P  xT  +  l)’1 
m  m  m  m 


I  *'  .  !  i.  >  A  T 

In  (9.1)  one  can  consider  the  term  (x^.  -  4_  x  )  x^  as  being  the 

m+l  m  m  m 

current  estimate  (based  on  the  transition  from  x  to  x  ..)  of  the 

m  m+i 

error  in  the  estimator;  4  .  Note  that  this  term  also  appears  in  (4.5). 

m 

Extending  this  line  of  reasoning,  the  factor  (1/Xm)  in  (9.1)  is  the- 

A 

weighting  given  this  error  term  in  forming  the  new  estimator,  4m+1 . 
Setting  X  *  1  will  Cause  every  error  so  generated  to  be  weighted 
equally  in  the  aggregate  estimate.  For  X<  1 ,  later  terms  will  be 
weighted  more  heavily  than  the  earlier  terms .  )•> 
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The  analogous  weighting  factor  in  (4.5)  is  the  term  a  P  .  The 

m  m 

two  formulations  are  seen  to  be  Identical  except  for  this  factor.  Thus, 
(9.1)  may  be  considered  to  be  a  special  case  of  (4.5),  with  non¬ 
optimum  weighting  of  the  error  term  (optimality  here  being  taken  in  the 
usual  least  squares  sense).  A  comparison  of  the  schemes  with  Xs  1 
is  presented  in  section  10. 

Lee  [  3  ]  considers  the  identification  of  a  system  with  a  scalar 
observable.  The  input  distribution  matrix  and  state  transition  matrix 
are  constrained  to  be  of  the  form 


I 


This  system  is  described  in  section  6.  Lee's  formulation  for  the  iden¬ 
tification  of  this  system  is  equivalent  to  that  herein  presented. 

t 

io.  Computational  results 

The  purpqse  of  this  section  is  to  obtain  computational  substan¬ 
tiation  6f  the  formulations  herein  presented.  To  this  end,  the  follow¬ 
ing  experimental  objectives  have  been  set  forth: 

1.  Test  die  algorithm  4  •  R,  R_1  for  convergence  to  *  . 

i  o 

2.  Test  the  recursive  algorithm  for  convergence. 

3.  Compare  l.  and  2.  above. 

4.  test  the  recursive  algorithm  with  additive  measurement  noise. 

5.  Investigate  convergence  under  varying  observability  con¬ 
straints  . 

6.  Compare  Ho's  [  2  ]  method  (equation  9.1)  with  the  recursive 
algorithm. 

Two  model  plants  were  selected  to  Implement  the  testing.  The 
first  is  a  simple  oscillator;  the  second  a  fourth  order  plant  with  two 
oscillatory  modes,  taken  to  be  representative  of  the  longitudinal 
dynamics  of  a  large  jet  transport  aircraft  during  normal  cruise.  [  1  ] 
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Parameters  for  the  two  system  are  as  follows: 

2n<*  order  plant:  w  =  1  rad/sec  T  ■  n/8  sec 

4^  order  plant:  ■  1.15  rad/sec  X ^  m  *35 

U>2  e  *11  rad/sec  $2  -  .035 

T  =  1.5  sec 

Testing  for  convergence .  (tests  1,  2  and  3) 

. . ‘  .  A 

All  the  experimentation  done  supported  the  contention  that  4*  does 
in  fact  converge  to  4  ,  and  that  the  estimator  4  is  therefore  unbiased. 
The  results  of  testing  with  the  2nc*  order  oscillator  are  presented  as 
being  typical.  In  fig.  (3)  ,  the  magnitude  of  the  elements  of  ^  are 
plotted  as  a  function  of  the  number  of  observations  used  in  forming  the 
estimator.  These  calculations  were  made  using  the  nonrecurslve  (batch 
processing)  algorithm. 

A 

In  fig.  (4)  ,  the  elements  of  4  are  plotted  as  a  function  of  the 

number  of  iterations,  and  were  calculated  using  the  recursive  algor- 

6 

ithm.  The  recursive  algorithm  was  initialized  by  setting  P.  ■  10  •  I. 

A  2  A 

Since  ft  was  demonstrated  in  section  7  that  cov(0.)  *  g.  P  ,  the 

l  l  m 

large  initial  value  for  demonstrates  this  uncertainty  in  the  initial 
value  of  4  1 ,  taken  to  be  the  identity  matrix.  The  resultant  estimator 
is  seen  to  be  very  nearly  equivalent  to  that  of  the  batch  processing 
technique,  by  a  comparison  of  figs.  (3)  and  (4).  Using  this  initial¬ 
ization  scheme  for  the  example  shown,  all  elements  of  the  recursive 
estimator  were  to  within  4  significant  figures  of  the  nonrecursive 
estimator  at  300  Iterations.  Identical  data  was  used  in  the  generation 
of  the  two  estimators . 

It  is  noted  that  the  recursive  form  of  the  algorithm  offers  several 
computational  advantages ,  including  its  suitability  for  real  time  im¬ 
plementation  and  the  fact  that  no  matrix  Inversion  is  required  in  its 
implementation.  Because  of  these  advantages,  the  recursive  form 
was  used  for  the  majority  of  the  remaining  testing ,  the  author  being 

A 

satisfied  that  no  degradation  of  the  estimator  4  would  result  from  so 
doing . 
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Testing  with  measurement  noise,  (test  4) 

For  the  remaining  tests,  especially  with  the  higher  order  system, 
it  is  desirable  to  introduce  a  scalar  measure  of  merit  for  the  estimators. 
While  the  definition  of  such  a  measure  is  certainly  arbitrary,  it  is  the 
author's  opinion  that  such  a  measure  should  make  a  comparison  of  the 
characteristic  roots  of  the  estimator  with  the  characteristic  roots  of 
the  actual  system.  This  method  seems  especially  appropriate  when 
comparing  matrices,  the  make  up  of  which  vary  considerably  with 
varying  observability  constraints.  By  definition,  then, 


where  the  X  are  the  characteristic  roots  of  the  matrix  4  ,  and  the 

*  *  A 

X 1  the  roots  of  the  estimator  4  . 

A  typical  result  of  estimating  with  uncorrelated  measurement 
noise,  using  the  4  order  model,  is  shown  in  Fig.  (5).  Here  for 
the  case  of  noisy  observations  is  compared  with  an  estimator  formed 
from  the  same  data,  but  without  the  additive  measurement  noise.  For 
this  example,  the  ratio  of  variances  for  excitation  and  measurement 
noise  was  taken  to  be 

4<a  *  /<*?  <  10  for  i,j  -  1,  2,  3,  4. 

wi  rJ 

Testing  with  constrained  observability,  (test  5.) 

A 

As  stated  in  section  8,  the  estimators  4*  ,  formed  from  systems 
having  constrained  observability,  in  general  are  not  unbiased.  How¬ 
ever,  several  interesting  properties  of  the  estimator  have  come  into 
evidence  from  computational  experimentation. 

The  first  experiment  under  constrained  observability  was  con- 
ducted  with  the  2  order  oscillator,  and  the  estimator  4*  formed 

from  observations  on  a  Single  mode  of  the  system . 

$.  •  .  ''  ' 
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Figure  7.  Comparison  of  convergence  for  estimators  formed  from  con¬ 
strained  and  unconstrained  observations.  2nd  order  system. 


Results  of  the  experiment  are  shown  in  figs .  (6)  and  (7) .  Convergence 
to  4*  ,  or  to  a  matrix  very  close  to  4  *  ,  is  demonstrated  in  this  case. 
Fig.  (7)  shows  a  comparison  of  the  constrained  observability  estimator. 

A  A 

4*  ,  and  the  estimator  formed  from  the  full  state  vector.  4  .  The 
measure  of  merit  suggested  by  Lee  [  3  ] 


M, 


*  -  in 

ii  *  ii  2 


,“6 


is  plotted  here.  At  300  interations.  M  equalled  10  w  for  the  case  of 

-4  £ 

the  vector  observable;  5x10  for  the  scalar  observable. 

Experimentation  with  the  4**1  order  system  under  constrained  ob¬ 
servability  failed  to  produce  a  convergent  estimator.  However,  five 
distinct  implementations  had  the  commonrproperty  of  correctly  iden¬ 
tifying  the  dominant  poles  of  the  system.  Fig.  (8)  shows  the  charac¬ 
teristic  roots,  at  300  iterations,  of  4  estimators,  based  respectively 
on  a  scalar,  2  vector,  3  vector,  and  full  state  vector  observable.  At 
2000  iterations,  the  roots  of  the  estimator  based  on  the  full  state 
vector  had  moved  to  within  3  significant  figures  of  the  actual  roots  of 
the  system ,  while  the  roots  of  the  other  three  estimators  had  not  moved 
appreciably  from  the  positions  plotted  in  fig.  (8). 

Fig.  (9)  compares  the  roots  of  an  estimator  formed  by  using  only 

th  V 

every  4  observation  with  an  estimator  formed/a  scalar  observable  at 

th 

each  sampling  instant.  Using  only  every  4  sample  is  the  method 
suggested  by  Lee  [  3  ]  ,  as  discussed  in  section  8.  Note  here  that 


the  state  transition  for  4  sampling  instants  is  given|4  *  (4T)  « 

4  4  2 

C  4  J  (T)  3  .  The  roots  of  4  J  (4T)  are  then  (  A  )  ,  where  are  the 

roots  of  4  *(T).  Making  this  calculation  places  the  dominant  poles  of 

- 

*(T)  in  close  proximity  of  the  poles  of  4,  as  shown  in  fig.  (10). 

*  S,  '  .  i  .  A 

In  section  8  it  was  suggested  that  the  order  of  the  estimator  4  ** 
for  the  case  of  constrained  observability  be  increased  to  2n-l  to  allow 

Oft','-  '  •••;  . '  '■  to  '  ,  ( 

for  the  inclusion  of  a  coloring  filter.  A  typical  result  of  such  an 

i  <  i 

augmented  estimator  is  shown  in  fig.  (11).  Also  shown,  for  comparison. 
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Z  PLANE  PLOT  OF  CHARACTERISTIC  ROOTS 


+  Roots  of  the  actual  state  transition  matrix 
o  Roots  of  ,  based  on  a  scalar  observable 

A  ^ 

^  Roots  of  4  ^  ,  based  on  a  2  vector  observable 
D  Roots  of  ^  *  ,  based  on  a  3  vector  observable 

A  ^ 

A  Roots  of  4  .  All  state  modes  observable 

•;v  ,  •  .  ,  ,,  ,  •  i.  .  *•  . 

Figure  8.  The  characteristic  roots  of  4  estimators  are  compared.  Only 

*  !  ’  V  J  ' 

the  estimator  4  is  converging  to  ♦.  The  dominant  (low  fre- 

’  *  :V  *0  ‘  *:>  •“’«—*  4  :  f  j  i  *  /'  V 

quency)  roots  of  all  estimators  are  within  3  significant 

A.-wlr  u  -ni  >  fif  »=:<{',  .  w- ;.f  -  .  »• 

figures  of  the  dominant  roots  of  ♦  .  300  iterations  were 

.  .’■■■■  .  >1/  ■  .  a.  '  :  u.  r 

used  to  generate  all  estimators . 
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Z  PLANE  PLOT  OF  CHARACTERISTIC  ROOTS 


+■  Roots  of  the  actual  state  transition  matrix 

O  Roots  of  ^  *  ,  based  on  a  scalar  observable 

A  Roots  of  4$  ,  based  on  a  scalar  observable, 

*  £|| 

using  only  every  4  observation 

Figure  9,  The  characteristic  roots  of  an  estimator  4  *  ,  formed  from 

th  1  A 

every  4  observation,  are  compared  with  the  roots  of  4  *, 

formed  from  successive  scalar  observations.  2000  iterations 

were  used  to  generate  the  estimators.  The  dominant  (low 

frequency)  roots  of  ^  ^  are  within  3  significant  figures  of 

the  roots  of  4  . 
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Z  PLANE  PLOT  OF  CHARACTERISTIC  ROOTS 


Figure  10.  The  characteristic  roots  of  4  *(T)  are  generated  from 
^£(4T).  These  roots  are  compared  with  the  roots  of 
The  data  of  fig.  (9)  Was  used  to  calculate  these 
roots . 


are  the  roots  of  a  4  th  order  estimator  based  on  a  scalar  observable. 

2000  Iterations  were  used  to  generate  the  data  shown. 

It  is  of  interest  to  note  the  proximity  of  the  roots  of  the  4**1  order 
estimator  to  4  of  the  roots  of  the  augmented  estimator. 

Comparison  with  Ho's  method,  (test  6.) 

Figure  (12)  shows  the  location  of  the  characteristic  roots  of  an 
estimator  generated  from  (9.1)  ,  the  formulation  proposed  by  Ho.  [2  ] 
These  roots  are  compared  with  the  roots  of  an  estimator  formed  from 
Identical  data,  using  (4.5),  the  recursive  algorithm.  300  iterations 
were  used  to  generate  these  estimators.  In  this  simulation,  the  roots 
of  the  estimator  generated  by  (9.1)  had  not  acquired  a  "convergent 
tendency"  at  300  iterations,  in  that  they  still  were  moving  quite  marked¬ 
ly  from  iteration  to  iteration.  In  contrast,  the  roots  of  the  estimator 
generated  by  (4 . 5)  were  confined  to  a  small  region  of  the  z  plane  after 
about  50  Iterations.  Thus,  it  is  not  intended  toimPly  that  the  estimator 
generated  from  Ho's  scheme  is  not  convergent,  but  simply  to  compare 
the  two  estimators  at  what  is  considered  a  typical  point  in  time.  Un¬ 
fortunately,  the  author  did  not  have  sufficient  time  to  accomplish  a 
more  rigorous  investigation  of  (9.1). 


Z  PLANE  PLOT  OF  CHARACTERISTIC  ROOTS 


i 


+  Roots  of  the  actual  state  transition  matrix 

A 

O  Roots  of  4  *  based  on  a  scalar  observable 
a  * 

X  Roots  of  ♦  **  ,  augmented  to  order  7 


Figure  11.  The  characteristic  roots  of  an  augmented  estimator 

A 

♦  **,  formed  from  a  scalar  observable,  are  compared 
with  the  roots  of  4  *.  2000  iterations  were  used  to 
generate  the  estimators . 
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Z  PLANE  PLOT  OF  CHARACTERISTIC  ROOTS 


•+*  Roots  of  the  actual  state  transition  matrix 
o  Roots  of  4  .  Full  state  vector  observable 
K  Roots  of  Ho's  estimator 


Figure  12.  The  characteristic  roots  of  an  estimator  proposed  by  Ho 

A 

are  compared  with  the  estimator  4  ,  formed  from  the  full 
state  vector.  300  iterations  were  used  to  generate  the 
estimators . 
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111-  ABSTRACT  * 

A  least  squares  estimator  % is  derived  for  the  state  transition  matrix  4- 

i 

of  a  linear,  stationary  sampled  data  system  operating  in  a  stochastic  environment. 

A 

The  estimator  A  is  shown  to  be  unbiased  and  minimum  variance  under  the  con¬ 
dition  of  full  observability  of  the  state  vector  of  the  system.  The  estimator  is 

<  ;  J 

also  Shown  to  be  the  Maximum  Likelihood  Estimator  for  the  case  of  the  stochastic 

i  l  i  f 

environment  having  Gaussian  statistics.  The  estimation  scheme  is  compared 

|  ,  !  I 

with  two  other  recently  published  estimation  schemes ,  both  of  which  are  shown 
to  be  special  cases  of  the  scheme  herein  presented. 
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